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c/3 A technique is described for constructing three-dimensional vector graphics 

i_i representations of planar regions bounded by cubic Bezier curves, such as 

_j smooth glyphs. It relies on a novel algorithm for compactly partitioning 

>■ planar Bezier regions into nondegenerate Coons patches. New optimizations 

^ are also described for Bezier inside-outside tests and the computation of 
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1. Introduction 
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Recent methods for lifting smooth two-dimensional (2D) font data into 

3 three dimensions (3D) have focused on rendering algorithms for the Graphics 

4 Processing Unit (GPU) [TJ]. However, scientific visualization often requires 

5 3D vector graphics descriptions of surfaces constructed from smooth font 

6 data. For example, while current CAD formats, such as the PDF-embeddable 

7 Product Representation Compact (PRC, precis in French) [2] format, allow 

8 one to embed text annotations, they do not allow text to be manipulated as 

9 a 3D entity. Moreover, annotations can only handle simple text; they are not 

10 suitable for publication-quality mathematical typesetting. 

11 In this work, we present a method for representing arbitrary planar re- 

12 gions, including text, as 3D surfaces. A significant advantage of this repre- 

13 sentation is consistency: text can then be rendered like any other 3D object. 

14 This gives one complete control over the typesetting process, such as kerning 

15 details, and the ability to manipulate text arbitrarily (e.g. by transformation 

16 or extrusion) in a compact resolution-independent vector form. In contrast, 

17 rendering and mesh-generation approaches destroy the smoothness of the 

18 original 2D font data. 

19 In focusing on the generation of 3D surfaces from 2D planar data, the 

20 emphasis of this work is not on 3D rendering but rather on the underly- 

21 ing procedures for generating vector descriptions of 3D geometrical objects. 

22 Vector descriptions are particularly important for online publishing, where 

23 no assumption can be made a priori about the resolution that will be used 

24 to display an image. As explained in Section [2| we focus on surfaces based 

25 on polynomial parametrizations rather than nonuniform rational B-splines 

26 (NURBS) P [19]. In Section |3] we describe a method for splitting an arbi- 

27 trary planar region bounded by one or more Bezier curves into nondegenerate 

28 Bezier patches. This algorithm relies on the optimized Bezier inside-outside 

29 test described in Section |4j The implementation of these algorithms in the 

30 vector graphics language Asymptote, along with the optimized 3D sizing 

31 algorithms presented in Section [5| is discussed in Section [6j 

32 Using a compact vector format instead of a large number of polygons to 

33 represent manifolds has the advantage of reduced data representation (essen- 

34 tial for the storage and transmission of 3D scenes) and the possibility, using 

35 relatively few control points, of exact or nearly exact geometrical descriptions 



36 of mathematical surfaces. For example, in Appendix Appendix A we show 



37 that a sphere can be represented to 0.05% accuracy with just eight cubic 



38 Bezier surface patches. 

39 2. Bezier vs. NURBS Parametrizations 

The atomic graphical objects in PostScript and PDF, Bezier curves and 
surfaces, are composed of piecewise cubic polynomial segments and tensor 
product patches, respectively. A segment 7(t) = Yli=o^i(^)^i ^^^ ^"-"^^ '^*-"^" 
trol points Pi, whereas a surface patch is defined by sixteen control points Pjj-: 

3 

(t{u,v) = {x{u,v),y{u,v)) = ^Bi{u)Bj{v)Pij. 

i,j=0 

40 Here Bi{u) = ( Ju*(l — m)'^~* is the ith cubic Bernstein polynomial. Just as a 

41 Bezier curve passes through its two end control points, a Bezier surface nec- 

42 essarily passes through its four corner control points. These special control 

43 points are called nodes. It is convenient to define the convex hull of a cubic 

44 Bezier segment or patch to be the convex hull (minimal enclosing polygon 

45 or polyhedron) of its control points. A straight segment is one in which the 

46 control points are colinear and the derivative of the Bezier parametrization 

47 is never zero (i.e. the control points are arranged in the same order as their 

48 indices) . 

49 It is often desirable to project a 3D scene to a 2D vector graphics for- 

50 mat understood by a web browser or high-end printer. Although NURBS 

51 are popular in computer-aided design [9] because of the additional degrees 

52 of freedom introduced by weights and general knot vectors, these benefits 

53 are tempered by both the lack of support for NURBS in popular 2D vector 

54 graphics formats (PostScript, PDF, SVG, EMF) and the algorithmic simpli- 

55 fications afforded by specializing to a Bezier parametrization. Bezier curves 

56 are also commonly used to describe glyph outlines. We therefore restrict 

57 our attention to (polynomial) Bezier curves and surfaces (even though both 

58 Asymptote and the 3D PRC format support NURBS). 

59 Unlike their Bezier counterparts, NURBS are invariant under perspective 

60 projection. This is only an issue if projection is done before the rendering 

61 stage, as is necessary when a 2D vector representation of a curve or surface 

62 is constructed solely from the 2D projection of its control points. It is there- 

63 fore somewhat ironic that NURBS are much less widely implemented in 2D 

64 vector graphics formats than in 3D. In 3D vector graphics applications, pro- 

65 jection to 2D is always deferred until rendering time, so that the invariance 



66 of NURBS under nonafiine projection is irrelevant. While NURBS provide 

67 exact parametrizations of familiar conic sections and quadric surfaces, non- 
68 trivial manifolds still need to be approximated as piecewise unions of under- 

69 lying exact primitives. We feel that the implement at ional simplicity of basic 

70 Bezier operations (computing subcurves and subsurfaces, points of tangency, 

71 normal vectors, bounding boxes, intersection points, arc lengths, and arc 

72 times) offsets for many practical applications the lower dimensionality of the 

73 Bezier subspace. 

74 3. Partitioning Curved 2D Regions 

75 In 3D graphics, text is often displayed with bit-mapped images, textures, 

76 or polygonal mesh approximations to smooth font character curves. To allow 

77 viewing of smooth text at arbitrary magnifications and locations, a nonpolyg- 

78 onal surface that preserves the curvature of the boundary curves is required. 

79 While it is easy to fill the outline of a smooth character in 2D, filling a 3D 

80 planar surface requires more sophisticated methods. One approach involves 

81 using surface filling algorithms for execution on CPUs [K]. When a vec- 

82 tor, rather than a rendered, image is desired, a preferable alternative is to 

83 represent the text as a parametrized surface. 

84 Methods based on common surface primitives in 3D modelling and ren- 

85 dering can be used to describe planar regions. One method trims the domain 

86 of a planar surface to the desired shape [T7]. While that approach is feasible, 

87 given adequate software support for trimming, this work describes a differ- 

88 ent approach, where each symbol is represented as a set of planar Bezier 

89 patches. We call this procedure bezulation since it involves a process similar 

90 to the triangulation of a polygon but uses cubic Bezier patches instead of 

91 triangles. To generate a surface representing the region bounded by a set of 

92 simple closed Bezier curves (intersecting only at the end points), algorithms 

93 were developed for (i) expressing a simply connected 2D region as a union of 

94 Bezier patches and (ii) breaking up a nonsimply connected region into sim- 

95 ply connected regions. (Selfintersecting curves can be handled by splitting at 

96 the intersection points.) These algorithms allow one to express text surfaces 

97 conveniently as Bezier patches. 

98 Bezulation of a simply connected planar region involves breaking the re- 

99 gion up into patches bounded by closed Bezier curves with four or fewer 

100 segments. This is performed by the routine bezulate (cf. Algorithm [1]) us- 

101 ing an adaptation of a naive triangulation algorithm, modified to handle 



curved edges, as illustrated in Figure [T} 



Input: simple closed curve C 
Output: array of closed curves A 
while C. segments > ^ do 
found ^ false; 
for n = 3 to 2 do 

for i = to C.segments-i do 

L ^ line segment between nodes i and i + n oi C; 

if countIntersections(^C,ivj = 2 and midpoint of L is 

inside C then 

p ^ subpath of C from node i to i + n; 
q ^ subpath of C from node i + n to i + C. segments; 
A.push(p+L); 
C ^ L + q; 
found ^ true; 
break; 
end 
end 

if found then 
I break; 
end 
end 
if not found then 

refine C by inserting an additional node at the parametric 
midpoint of each segment; 
end 
end 
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Algorithm 1: bezulate partitions a simply connected region. 

A line segment lies within a closed curve when it intersects the curve 
only at its endpoints and its midpoint lies strictly inside the curve. If after 
checking all connecting line segments between nodes separated by n = 3 or 
n = 2 segments, none of them lie entirely inside the shape, the original curve 
is refined by dividing each segment of the curve at its parametric midpoint. 
The bezulation process then continues with the refined curve. This algorithm 
can be modified to subdivide more optimally, for example, to avoid elongated 
patches that sometimes lead to rendering problems. 

If the region is convex. Algorithm [l] is easily seen to terminate: all con- 



(a) 



i + l 




(b) 



(c) 



(d) 



(e) 



Figure 1: The bezulate algorithm. Starting with the original curve (a), several possible 
connecting line segments (shown in red) between nodes separated by n = 3 or n = 2 
segments are tested. Connecting line segments are rejected if they do not lie entirely 
inside the original curve. This occurs when the midpoint is not inside the curve (b) or 
when the connecting line segment intersects the curve more than twice (c) . If a connecting 
line segment passes both tests, the shaded section is separated (d) and the algorithm 
continues with the remaining curve (e). 
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necting line segments are admissible, and each patch removal decreases the 
number of points in the curve. Moreover, from the point of view of Algo- 
rithm [T| upon sufficient subdivision a non-convex region eventually becomes 
indistinguishable from a polygon, in which case the algorithm reduces to a 
straightforward polygonal triangulation. 

3.1. Nonsimply Connected Regions 

Since the bezulate algorithm requires simply connected regions, nonsim- 
ply connected regions must be handled specially. The "holes" in a nonsimply 
connected domain can be removed by partitioning the domain into a set of 
simply connected regions, each of which can then be bezulated. 

For convenience we define a top-level curve to be a curve that is not 
contained inside any other curve and an outer (inner) curve to be the outer 
(inner) boundary of a filled region. With these definitions, the glyph "%" 
has two inner curves and two top-level curves that are also outer curves. 

The algorithm proceeds as follows. First, to determine the topology of 
the region, the curves are sorted according to their relative insidedness, as 
determined by the nonzero winding number rule. Since the curves are as- 
sumed to be simple, any point on an inner curve can be used to test whether 
that curve is inside another curve. The result of this sorting is a collection 
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of top-level curves grouped with the curves they surround. Each of these 
groups is treated independently. 

Figure [3] illustrates the partition routine (cf. Algorithm [2]) . Each group 
is examined recursively to identify regions bounded by inner and outer curves. 
First, the inner curves in the group are sorted topologically to find the inner 
curves that are top-level curves with respect to the other inner curves. The 
inner curves that are not top-level curves are processed with a recursive call to 
partition. The nonsimply connected region between the outer (top-level) 
curve and the inner (top-level) curves is now split into simply connected 
regions. This is illustrated in Figure |2} The intersections of the inner and 
outer curves with a line segment from a point on an inner curve to a point 
on the outer curve are found (either via subdivision or a numerically robust 
cubic root solver). Consecutive intersections of this line segment, at points 
A and B, on the inner and outer curves, respectively, are selected. Let t^ 
be the value of the parameter used to parameterize the outer curve at B. 
Starting with A = 1, A is halved until the line segment AC, where C is 
the point on the outer curve at t^ + A, does not intersect the outer curve 
more than once, does not intersect any inner curve (other than once at A), 



150 and the region bounded by AB, AC, and BC does not contain any inner 

151 curves. Once A and the point C have been found, the outer curve, less 

152 the segment between B and C, is merged with BA, followed by the inner 



curve and then AC. The region bounded by AB, AC, and BC is a simply 
connected region. Additional simply connected regions are found when the 
outer curve is merged with the other inner curves. Once the merging with 
all inner curves has been completed, the outer curve becomes the boundary 
of the final simply connected region. 

The recursive algorithm for partitioning nonsimply connected regions into 
simply connected regions is summarized below. The function sort returns 
groups of top-level curves and the curves they contain. However, it is not 
recursive; the inner curves are not sorted. The function merge returns the 
simply connected regions formed from the single outer curve and multiple 



(a) 



(b) 



(c) 



(d) 



(e) 



(f) 



Figure 2: Splitting of non-simply connected regions into simply connected regions. Starting 
with a non-simply connected region (a), the intersections between each curve and an 
arbitrary line segment from a point on an inner curve to the outer curve are found (b). 
Consecutive intersections of this line segment, at points A and B, on the inner and outer 
curves, respectively, identify a convenient location for extracting a region. One searches 
along the outer curve for a point C such that the line segment AC intersects the outer 
curve no more than once, intersects an inner curve only at A, and determines a region 
ABC between the inner and outer curves that does not contain an inner curve. Once such 
a region is found (c), it is extracted (d). This extraction merges the inner curve with the 
outer curve. The process is repeated until all inner curves have been merged with the 
outer curve, leaving a simply connected region (e) that can be split into Bezier surface 
patches. The resulting patches and extracted regions are shaded in (f). 



163 inner curves that are supplied to it. 



Input: array of simple closed curves C 

Output: array of closed curves A 

foreach group of nested curves G in sort(^Cj do 

innerG roups ■(— sort(G. inner Curves); 

foreach group of nested curves H in innerGroups do 
I A.push(partition(H.innerCurves)); 

end 

A.push(merge(G.toplevel, top-level curves of all groups in 

innerGroups)); 
end 
return A; 

Algorithm 2: partition splits nonsimply connected regions into sim- 
ply connected regions. The pseudo-code functions sort and merge are 
described in the text. 
The routines bezulate and partition were used to typeset the T^^X 




merge, bezulate 



merge, bezulate 



bezulate 




Figure 3: Illustration of the partition algorithm. The five curves that define the outlines 
of the Greek characters a and are passed in a single array to partition. 




Figure 5: Zoomed view of Figure |4] gener- 
Figure 4: Apphcation of the bezulate and ated from the same vector graphics data, 
partition algorithms to hft the Gaussian The smooth boundaries of the characters 
integral to three dimensions. emphasize the advantage of a 3D vector font 

description. 






Figure 6: Subpatch boundaries for Figurel4]as determined by the bezulate and partition 
algorithms. 
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165 in the interactive 3D diagram shown in Figure |4] and magnified, to emphasize 

166 the smooth font boundaries, in Figure [5j The computed subpatch boundaries 

167 are indicated in Figure [6j 



168 Figure A.12 in Appendix Appendix A illustrates how bezulate is used 



169 in mathematical drawings to lift T[t;X to three dimensions. Referring to 
the interactive 3D PDF version of this articlq3 one see that the labels in 



171 Figure A.12 have been programmed to rotate interactively so that they always 

172 face the camera; this feature, implemented with Javascript, is known as 

173 billboard interaction. 

174 Developing Bezier versions of more sophisticated triangulation algorithms 

175 would be an interesting future research project. The rendering technique 

176 of Ref. [IS] could be modified to produce Bezier patches, but this would 

177 produce more patches than bezulate. For example, the "e" shown in Fig. 3 

178 of Ref. [15] corresponds to roughly twice as many (4-segment) patches as 

179 the ten patches generated by bezulate for the "e" in Fig. [6] Since our 
ISO interest is in compact 3D vector representations, the objective of this work 
181 is to minimize the number of generated patches. In contrast, in real-time 

132 rendering, one aims to minimize the overall execution time. 

133 3.2. Nondegenerate Planar Bezier Patches 

134 The bezulate algorithm described previously decomposes regions bounded 

135 by closed curves (according to the nonzero winding number rule) into subre- 

136 gions bounded by closed curves with four or fewer segments. Further steps are 

137 required to turn these subregions into nondegenerate Bezier patches. First, if 

138 the interior angle between the incoming and outgoing tangent directions at a 

139 node is greater than 180°, the boundary curve is split at this node by follow- 

190 ing the interior angle bisector to the first intersection with the path. This is 

191 done to guarantee that the patch normal vectors at the nodes all point in the 

192 same direction. Next, curves with less than four segments are supplemented 

193 with null segments (four identical control points) to bring their total number 

194 of segments up to four. A closed curve with four segments defines the twelve 



^See http : //asymptote . sourceforge .net/axticles/ 
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195 boundary control points of a Bezier patch in the x-y plane. The remaining 

196 four interior control points {-Pn, -P12, -P21, -P22} are then chosen to satisfy the 

197 Coons interpolation [3 [TOl QJ 

3 

(t{u, v) = Y, [(1 - v)Bi{u)P,fl + vB,{u)Pi^^ + (1 - u)Bi{v)Po,i + uBi{v)P^^i] 

-(1 - n)(l - t;)Po,o - (1 - w)^Po,3 - «(1 - w)P3,o - «t^P3,3. 

The resulting mapping (t{u,v) need not be bijective [22l [Ml [H], even if 
the corner control points form a convex quadrilateral (despite the fact that 
a Coons patch for a convex polygon is always nondegenerate). In terms of 
the 2D scalar cross product pxq = pxQy — Vylx, the Coons patch is seen to 
be a diffeomorphism of the unit square D = [0, 1] x [0, 1] if and only if the 
Jacobian 

J{u,v) = ^jf^ = VuXXV,y= J2 Bl{u)B,{v)B,{u)B',{v)PijXPki 

i,j,k,i=0 

198 (the z component of the corresponding 3D normal vector) is sign-definite. 

199 Since J{u, v) is a continuous function of its arguments, this means that J 

200 must not vanish anywhere on D. A sign reversal of the Jacobian can manifest 

201 itself as an outright overlap of the region bounded by the curve or as an 

202 internal multivalued wrinkle, as illustrated in Figure [7| Rendering problems, 

203 such as the black smudges visible in Figures [TJ^b) and (e), can occur where 

204 isolines collide. 

Randrianarivony and Brunnett [22] (and later H. Lin et al. [H]) describe 
sufficient conditions for J{u,v) to be nonzero throughout D. In the case of 
a cubic Bezier patch, the 36 quantities 

i+k=p j+e=q \/\/\J/\/ 

205 where t/jj = Pi+ij — Pij and Vij = Pij+i — Pij, are required to be of 

206 the same sign. This follows from the fact that J{u,v) = '^pq^oTpgU^v'^^l — 

207 uf-P{l-vf-'^. 

208 Randrianarivony et al. show further that every degenerate Coons patch 

209 can be decomposed into a finite union of nondegenerate subpatches (some 

210 with reversed orientation). However, the adaptive subdivision algorithm they 

12 







(d) 



(e) 



(f) 



Figure 7: Degeneracy in a Coons patch. The dots indicate corner control points (nodes) 
and the open circles indicate the points of greatest degeneracy on the boundary, as deter- 
mined by the quartic root solver: (a) overlapping isoline mesh; (b) overlapping patch; (c) 
nonoverlapping subpatches; (d) internally degenerate isoline mesh; (e) internally degener- 
ate patch; (f) nondegenerate subpatches. 
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211 propose to exploit this fact does not prescribe an optimal boundary point 

212 at which to do the splitting. A better algorithm is based on the following 

213 elementary theorem, which provides a practical means of detecting Coons 

214 patches with degenerate boundaries. 

Theorem 1 (Nondegenerate Boundary). Consider a closed counter-clockwise 
oriented four- segment curve p in the x^ plane such that the interior angles 
formed by the incoming and outgoing tangent vectors at each node are less 
than or equal to 180°. Let J{u, v) be the Jacobian of the corresponding Coons 
patch constructed from p, with control points Pij, and define the fifth-degree 
polynomial 

3 

i,j=0 

215 If f{u) > whenever f'{u) = on u E (0,1), then J{u,0) > on [0,1]. 

216 Otherwise, the minimum value of J{u,0) occurs at a point where f'{u) = 0. 

Proof First we note, since B[{0) = -fio(O) = 3 and 5^(0) = S^(0) = 0, 
that J{u, 0) = 3/(m) and 

J(0, 0) = 3/(0) = 9(Pi,o - Po,o) X (Po,i - Pofi) > 

217 since this is the cross product of the outgoing tangent vectors at Po,o- Like- 

218 wise, J(l, 0) = 3/(1) > 0. We know that the continuous function / must 

219 achieve its minimum value on [0,1] at some -u G [0, 1]. If / were negative 

220 somewhere in (0, 1) we could conclude that f{u) < 0, so that u G (0, 1), and 

221 hence / would have an interior local minimum at u, with f'{u) = 0. But this 

222 is a contradiction, given that f{u)>0 whenever f'{u) = 0. D 

The significance of Theorem [T] is that it affords a means of detecting a 
point u on the boundary where the Jacobian is most negative. This requires 
finding roots of the quartic polynomial 



fiu) = [B':{u)B,{u) + Bi{u)B'{u)]P,,oX{P, 



j,i - JTjfi) 



The coefficients of this quartic polynomial can be computed using the polyno- 
mials Mij = [B'l B j + B[B'j) / ?) tabulated in Table fll The method of Neumark 
|16j . which relies on numerically robust cubic and quadratic root solvers, is 
then used to find algebraically all real roots of the quartic equation f'{u) = 
that lie in (0, 1). The Jacobian is computed at each of these points; if it is 
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/ 5 - 20u + 30«2 - 20«3 + 5m'' -3 + 24u - ^Au'^ + A8u^ - I5u'^ -6u + 27u^ ~ 36u^ + 15u^ -3u'^ + 8u^ - 5u'^ \ 

-7 + 36«-66w2 + 52u3 _ i5„4 3 _ 36m + lOSw^ _ 120m^ + 45m'' 6u - Abu"^ + 84u^ - 45u'^ Sm^ - IBm^ + 15m4 

2 - 18m + 45m2 - 44m3 + ISm* 12m - 63m2 + 96u^ - 45m'' ISm^ - 60u^ + 45m4 Sm^ - 15m"' 

V 2m - 9m2 + I2u3 _ 5^4 g„2 _ 24„3 ^ ^^5^4 i2u^ - ISm" 5m'' J 

Table 1: Coefficients of the polynomials M^ = {Bl'Bj + B[B'A/i. 



228 negative anywhere, the point where it is most negative is determined. The 

229 patch is then spht along an interior hne segment perpendicular to the tan- 

230 gent vector at this point. The next intersection point of the patch boundary 

231 with this line is used to split the patch into two pieces. Each of these pieces 

232 is then treated recursively (beginning with an additional call to bezulate, 

233 should the new boundary curve happen to have five segments). 

234 If a patch possesses only internal degeneracies, like the one in Figure [7](d), 

235 the patch boundary is arbitrarily split into two closed curves, say along the 

236 perpendicular to the midpoint of some nonstraight side. The blue lines in 

237 Figures [TJ^b) and (f ) illustrate such a midpoint splitting. The arguments of 

238 Randrianarivony et al. [22j establish that only a finite number of such sub- 

239 divisions will be required to obtain a nondegenerate patch. Nondegenerate 

240 subpatches oriented in the direction opposite to the normal vector corre- 

241 sponding to the original oriented curve should be discarded to avoid rendering 

242 interference with correctly aligned overlying subpatches. 

243 The blue lines in Figure [TFc) show that our quartic algorithm generates 

244 six subpatches, a substantial improvement over the nine subpatches produced 

245 by adaptive midpoint subdivision [22] in Figure [TJ^b). Figure [TJ^c) also em- 

246 phasizes the ability of the quartic root algorithm to detect the optimal (most 

247 degenerate) points (circled) for splitting the boundary curve. As mentioned 

248 earlier, in both cases, it is possible that splitting can lead to curves with five 

249 segments. Such curves are split further by the bezulate algorithm so that 

250 any degeneracy of the resulting subpatches can be addressed. 

251 Since an algebraic quartic root solver is an explicit algorithm, optimal 

252 subdivision of patches introduces minimal overhead compared to adaptive 

253 midpoint subdivision. In our implementation, the costs of adaptive mid- 
254 point subdivision for Figures [Tf^b) and Figure [7](f ) were approximately the 
255 same. Using optimal subdivision in Figure [7](c) was 34% faster than adaptive 
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256 midpoint splitting, whereas there was only 2% additional overhead in check- 

257 ing for boundary degeneracy in Figure [71(f) (which possesses only internal 

258 degeneracy). Patches having only internal degeneracy arise relatively rarely 

259 in practice, but when they do, the subpatches obtained by adaptive midpoint 

260 subdivision also tend to exhibit internal degeneracy. Once internal degener- 

261 acy has been detected in a patch, we find that it is typically more efficient 

262 not to check its degenerate subpatches for boundary degeneracy (otherwise 

263 the overhead in checking for boundary degeneracy in Figure [Wf ) would grow 

264 to 50%). Of course, since our interest is not in real-time rendering but in 

265 surface generation, the real advantage of optimal subdivision is that it can 

266 significantly reduce the number of generated patches (e.g. Figure JTJ^c) has 
one-third fewer patches than Figure mb)). 



267 



268 4. An Optimized Bezier Inside— Outside Test 

269 Although PostScript has an infill function for testing whether a par- 

270 ticular point would be painted by the PostScript fill command, this is only 

271 an approximate digitized test corresponding to the resolution of the output 

272 device. Our bezulate routine requires a vector graphics algorithm, one that 

273 yields the winding number of an arbitrary closed piecewise Bezier curve about 

274 a given point. 

275 A straightforward generalization of the standard ray-to-infinity method 

276 for computing winding numbers of a polygon about a point requires the so- 

277 lution of a cubic equation. As is well known, the latter problem can become 

278 numerically unstable as two or three roots begin to coalesce. While a con- 

279 ventional ray-curve (or ray-patch) intersection algorithm based on recursive 

280 subdivision [T7j could be employed to count intersections by actually finding 

281 them, this typically entails excessive subdivision. 

282 A more efficient but still robust subdivision method for computing the 

283 winding number of a closed Bezier curve arises from the topological obser- 

284 vation that if a point z lies outside the convex hull of a Bezier segment, the 

285 segment can be continuously deformed to a straight line segment between its 

286 endpoints, without changing its orientation relative to the point z. A given 

287 point will typically lie outside the convex hull of most segments of a Bezier 

288 curve. The orientation of these segments relative to the given point can be 

289 quickly and robustly determined, just as in the usual ray method for poly- 

290 gons, to determine the contribution, if any, to the winding number. For this 
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Figure 8: The BezierWindingNumber algorithm. Since z Ues inside the convex huU of one 
Bezier segment, indicated by the Ught shaded region, that segment must be subdivided. 
On subdivision, z now Hes outside the convex huUs of the subsegments, indicated by the 
dark shaded regions; these subsegments may be continuously deformed to straight line 
segments between their endpoints, without crossing z. The usual polygon inside-outside 
test may then be applied: the green ray establishes a winding number contribution of +1 
due to the orientation of z with respect to the blue line. 



291 purpose, Jonathan Shewchuk's public-domain adaptive precision predicates 

292 for computational geometry [23] are highly recommended. 

293 In the infrequent case where z lies on or inside the convex hull of a seg- 

294 ment, de Casteljau subdivision is used to split the Bezier segment about 

295 its parametric midpoint. Typically the convex hulls of the resulting sub- 

296 segments will overlap only at their common control point, so that z can lie 

297 strictly inside at most one of these hulls. This observation is responsible 

298 for the efficiency of the algorithm: one continues subdividing until the point 

299 is outside the convex hull of both segments or until machine precision is 

300 reached, as illustrated in Figure |8| 

301 The orientation of segments whose convex hulls do not contain z can be 

302 handled by using the topological deformation property together with adap- 

303 tive precision predicates. Denoting by straightContribution(P,Q,z) the 

304 usual ray method for determining the winding number contribution of a line 

305 segment PQ relative to a point z, the contribution from a Bezier segment S 
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can be computed as curvedContribution(S,z) (Algorithm 3). 


Input: segment S, pair z 


Output: winding number contribution of S about z 


W^O; 


if z lies within or on the convex hull of S then 




foreach subsegment s of S do 




W ^ W + curvedContribution(s,z); 




end 


else 


W ^ W + straightContribution(S.beginpoint,S.endpoint,z); 


end 


return W; 



Algorithm 3: curvedContribution(S,z) determines the winding 
number contribution from a Bezier segment S about z. 
The winding number for a closed curve p about z may then be evaluated 
with the algorithm bezierWindingNumber(C,z) (Algorithm El) . 



Input: curve C, pair z 

Output: winding number of C about z 

foreach segment S of C do 
if S is straight then 

I W ^— W + straightContribution(S.beginpoint,S.endpoint,z) 
else 

I W ^ W + curvedContribution(S,z); 
end 
end 
return W; 



Algorithm 4: bezierWindingNuinber(C,z) computes the winding 

number of a closed Bezier curve C about z. 
A practical simplification of the above algorithm is the widely used op- 
timization of testing whether a point is inside the 2D bounding box of the 
control points rather than their convex hull. Since the convex hull of a Bezier 
segment is contained within the bounding box of its control points, one can 
replace "convex hull" by "control point bounding box" in the above algorithm 
without modifying its correctness. One can easily check numerically that the 
cost of the additional spurious subdivisons is well offset by the computational 



318 savings in testing against the control point bounding box. 

319 5. Global Bounds of Directionally Monotonic Functions 

320 We now present efficient algorithms for computing global bounds of real- 

321 valued directionally monotonic functions / : M'^ — > M defined over a Bezier 

322 surface cr{u,v). By directionally monotonic we mean that the restriction 

323 of / to each of the three Cartesian directions is a monotonic function; if / 

324 is differentiable this means that / has sign-semidefinite partial derivatives. 

325 These algorithms can be used to compute the 3D bounding box of a Bezier 

326 surface, the bounding box of its 2D projection, or the optimal field-of-view 

327 angle for sizing a 3D scene (cf. Fig. [o]). The key observation is that the 

328 convex hull property of a Bezier patch holds independently in each direction 

329 and even under inversions like z — )■ 1/z. 

330 A naive approach to computing the bounding box of a Bezier patch re- 

331 quires subdivision whenever the 3D bounding boxes overlap in any of the 

332 three Cartesian directions. However, the number of required subdivisons 

333 can be greatly reduced by decoupling the three directions: in Algorithm |5} 

334 the problem is split into finding the maximum and minimum of the three 

335 Cartesian axis projections f{x,y,z) = x, f{x,y,z) = y, and f{x,y,z) = z 

336 evaluated over the patch. This requires a total of six applications of Al- 

337 gorithm [5j By convexity, the extrema of these special choices for / over a 

338 convex polyhedron C occur at vertices of C. 

339 More general choices of directionally monotonic functions / are also of 

340 interest. For example, to determine the bounding box of the 2D perspective 

341 projection (based on similar triangles) of a surface, one can apply Algorithmjo] 

342 in eye coordinates to the functions f{x, y, z) = x/z and f{x, y, z) = y/z. This 

343 is useful for sizing a 3D object in terms of its 2D projection. For example, 

344 these functions were used to calculate the optimal field-of-view angle 13.4° 

345 for the Klein bottle shown in Figure [9j 

346 For an arbitrary directionally monotonic function /, we note that 

acC^ f{a) C /(C). (1) 

347 Our algorithms exploit Eq. ([I| together with de Casteljau's subdivision 

348 algorithm and the fact that a Bezier patch is confined to the convex hull of 

349 its control points. However, a patch is only guaranteed to intersect its convex 

350 hull at the four corner nodes. 
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Figure 9: A Bezier approximation to a projection of a four-dimensional Klein bottle to 
three dimensions. The FunctionMax algorithm was used to determine the optimal field 
of view for this symmetric perspective projection of the scene from the camera location 
(25.09,-30.33,19.37) looking at (-0.59,0.69,-0.63). The extruded 3D TgX equations 
embedded onto the surface provide a parametrization for the surface over the domain 
uxv e [0,27r] X [0,27r]. 
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351 For the special case where / is a projection onto the Cartesian axes, the 

352 function CartesianMax(/, P, /(Pqo); c^) given in Algorithm IS] computes the 

353 global maximum M of a Cartesian axis projection / : M^ — )■ M over a Bezier 

354 patch P to recursion depth d. Here, the value /(-Poo) provides a convenient 

355 starting value (lower bound) for M; if the maximum of a surface consisting 

356 of several patches is desired, the value of M from previous patches is used 

357 to seed the calculation for the subsequent one. The algorithm exploits the 

358 fact that the extrema of each coordinate over the convex hull C of P occur 

359 at vertices of C. First, one replaces M by the maximum of / evaluated at 

360 the four corner nodes and the previous value of M. If the maximum of the 

361 function evaluated at the remaining 12 control points is less than or equal 

362 to M, the subpatch can be discarded (by Eq. [l| noting that the maximum 

363 of f{C) occurs at a control point and hence cannot exceed M). Otherwise, 

364 the patch is subdivided along the u = v = 1/2 isolines and the process is 

365 repeated using the new value of M. The method quickly converges to the 

366 global maximum of / over the entire patch. 



Input: real function f(triple), patch P, 
Output: real M 

M ^ max(M, f (Pqo), f (Pos), f (i'so), f (i^ 
if depth = then 

return M; 
end 


real M, integer 
•33)); 


depth 


V ^ max(f(Poi),f(Po2),f(Pio),f(Pii) 

f(P2o),f(P2l),f(P22),f(P23) 

if V < M then 


f(Pl2),f(Pl3), 
f(P3l),f(P32)) 




return M; 
end 








foreach subpatch S of P do 

M ^ max(M, FunctionMax(f, S, M, d 
end 


epth - 


-1)); 




return M; 









Algorithm 5: CartesianMax(f,P, M , depth) returns the maximum of 
M and the global bound of a Cartesian component f of a Bezier patch 
P evaluated to recursion level depth. 

367 For a general directionally monotonic function / (consider f{x, y, z) = xy 
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368 over C = d{{x,y,0) :0<x<l, < y < x}), the maximum of /(C) need 

369 not occur at vertices of C: one instead needs to examine the function value 

370 at the appropriate vertex of the bounding box of C. For example, if / is a 

371 monotonic increasing function in each of the three Cartesian directions, 

Ccbox(a,b)^/(C)c[/(a),/(b)], (2) 

372 where box(a, b) denotes the 3D box with minimal and maximal vertices a 

373 and b, respectively. 

374 The global maximum M of a directionally monotonic increasing function 

375 / : M^ — i- M over a Bezier patch P can then be efficiently computed to 

376 recursion depth d by calling the function FunctionMax(/, P, /(Pqo); d) given 

377 in Algorithm |6| First, one replaces M by the maximum of / evaluated at 

378 the four corner nodes and the previous value of M. One then computes the 

379 vertex b of the bounding box of the convex hull C of P. If the maximum of 

380 the function evaluated at b is less than or equal to M, the subpatch can be 

381 discarded. Otherwise, the patch is subdivided along the u = v = 1/2 isolines 

382 and the process is repeated using the new value of M. 

383 6. 3D Vector Typography 

384 Donald Knuth's T[t;X system [13], the de-facto standard for typesetting 

385 mathematics, uses Bezier curves to represent 2D characters. TgX provides 

386 a portable interface that yields consistent, publication quality typesetting of 

387 equations, using subtle spacing rules derived from centuries of professional 

388 mathematical typographical experience. However, while it is often desirable 

389 to illustrate abstract mathematical concepts in TgX documents, no compati- 

390 ble descriptive standard for technical mathematical drawing has yet emerged. 

391 The recently developed Asymptote languagqj aims to fill this gap by 

392 providing a portable TgX-aware tool for producing 2D and 3D vector graph- 

393 ics [5]. In mathematical applications, it is important to typeset labels and 

394 equations with TgX for overall consistency between the text and graphical el- 

395 ements of a document. In addition to providing access to the TgX typesetting 

396 system in a 3D context, ASYMPTOTE also fills in a gap for nonmathematical 



^available from http://asymptote.sourceforge.net under the GNU Lesser General 
Public License. 
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Input: real function f(triple), patch P, real M, integer depth 
Output: real M 

M^max(M,f(Poo),f(Po3),f(P3o),f(P33)); 
if depth = then 

I return M; 
end 

X ^ max(A ■ Pij '■ <i,j < 3); 
y ^ max(^ ■ P^ '■ < i,j < 3); 
z ^ max(£ • P^ '■ < i,j < 3); 
if /((x,y,z)) < M then 

I return M; 
end 
foreach subpatch S of P do 

I M ^ max(M, FunctionMax(f, S, M, depth - 1)); 
end 
return M; 

Algorithm 6: FunctionMax(f,P,M, depth) returns the maximum of M 
and the global bound of a real- valued directionally monotonic increasing 
function f over a Bezier patch P evaluated to recursion level depth. Here 
X, y, z are the Cartesian unit vectors. 
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397 applications. While open source 3D bit-mapped text fonts are widely avail- 

398 ablejj resources currently available for scalable (vector) fonts appear to be 

399 quite limited in three dimensions. 

400 Asymptote was inspired by John Hobby's METRPOST (a modified ver- 

401 sion of METRFONT, the program that Knuth wrote to generate the TgX 

402 fonts), but is more powerful, has a cleaner syntax, and uses IEEE floating 

403 point numerics. An important feature of Asymptote is its use of the simplex 

404 linear programming method to solve overall size constraint inequalities be- 

405 tween fixed-sized objects (labels, dots, and arrowheads) and scalable objects 

406 (curves and surfaces). This means that the user does not have to scale man- 
ually the various components of a figure by trial-and-error. The 3D versions 
of Asymptote's deferred drawing routines rely on the efficient algorithms 
for computing the bounding box of a Bezier surface, along with the bounding 
box of its 2D projection, described in Sec. [5j Asymptote natively generates 

411 PostScript, PDF, SVG, and PRC |2] vector graphics output. The latter is a 

412 highly compressed 3D format that is typically embedded within a PDF file 

413 and viewed with the widely available Adobe Reader software. 

414 The biggest obstacle that was encountered in generalizing Asymptote 

415 to produce 3D interactive output was the fact that TgX is fundamentally a 

416 2D program. In this work, we have developed a technique for embedding 

417 2D vector descriptions, like TgX fonts, as 3D surfaces (2D vector graphics 

418 representations of TgK output can be extracted with a technique like that 

419 described in Ref. [6]). While the general problem of filling an arbitrary 3D 

420 closed curve is ill-posed, there is no ambiguity in the important special case 

421 of filling a planar curve with a planar surface. 

422 Since our procedure transforms text into Bezier patches, which are the 

423 surface primitives used in Asymptote, all of the existing 3D Asymptote 

424 algorithms can be used without modification. Together with the 3D gen- 

425 eralization of the METRFONT curve operators described by |H |5], these 

426 algorithms comprise the 3D foundation for the Tj^jX-aware vector graphics 

427 language Asymptote. 

428 6.1. 3D Arrowheads 

429 Arrows are frequently used in illustrations to draw attention to important 

430 features. We designed curved 3D arrowheads that can be viewed from a 



For example, see http://www.opengl.org/resources/features/fontsurvey/ 
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wide range of angles. For example, the default 3D arrowhead was formed by 
bending the control points of a cone around the tip of a Bezier curve. Planar 
arrowheads derived from 2D arrowhead styles are also implemented; they are 
oriented by default on a plane perpendicular to the initial viewing direction. 

The bezulate 



435 Examples of these arrows are displayed in Figures 10 and 11 



algorithm was used to construct the upper and lower faces of the filled (red) 



planar arrowhead in Fig. 11 




Figure 10: Three-dimensional revolved arrowheads in Asymptote. 




Figure 11: Planar arrowheads in Asymptote. 



439 
440 
441 
442 
443 
444 
445 
446 
447 
448 
449 



6.2. Double Deferred Drawing 

Journal size constraints typically dictate the final width and height, in 
PostScript coordinates, of a 2D or projected 3D figure. However, it is often 
convenient for users to work in more physically meaningful coordinates. This 
requires deferred drawing: a graphical object cannot be drawn until the actual 
scaling of the user coordinates (in terms of PostScript coordinates) is known 
[5j. One therefore needs to queue a function that can draw the scaled object 
later, when this scaling is known. Asymptote's high-order functions provide 
a flexible mechanism that allows the user to specify either or both of the 3D 
model dimensions and the final projected 2D size. This requires two levels of 
deferred drawing, one that first sizes the 3D model and one that scales the 
resulting picture to fit the requested 2D size [6]. The 3D bounding box of 
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450 a Bezier surface, along with the bounding box of its 2D projection, can be 

451 efficiently computed with the method described in Section [5] 

452 6.3. Efficient Rendering 

453 Efficient algorithms for determining the bounding box of a Bezier patch 

454 also have an important application in rendering. Knowing the bounding box 

455 of a Bezier patch allows one to determine, at a high level, whether it is in the 

456 field of view: offscreen Bezier patches can be dropped before mesh generation 
occurs |ljy . This is particularly important for a spatially adaptive algorithm 
as used in Asymptote's OPENGL-based renderer, which resolves the patch 

459 to one pixel precision at all zoom levels. Moreover, to avoid subdivision 

460 cracks, renderers typically resolve visible surfaces to a uniform resolution. 

461 It is therefore important that offscreen patches do not force an overly fine 

462 mesh within the viewport. As a result of these optimizations, the native 

463 Asymptote adaptive renderer is typically comparable in speed with the 

464 fixed-mesh PRC renderer in Adobe Reader, even though the former yields 

465 higher quality, true vector graphics output. 

466 7. Conclusions 

467 In this work we have developed methods that can be used to lift smooth 

468 fonts, such as those produced by T[t;X, into 3D. Treating 3D fonts as sur- 

469 faces allows for arbitrary 3D text manipulation, as illustrated in Figures [5] 

470 and[9J The bezulate algorithm allows one to construct planar Bezier surface 

471 patches by decomposing (possibly nonsimply connected) regions bounded by 

472 simple closed curves into subregions bounded by closed curves with four or 

473 fewer segments. The method relies on an optimized subdivision algorithm 

474 for testing whether a point lies inside a closed Bezier curve, based on the 

475 topological deformation of the curve to a polygon. We have also shown how 

476 degenerate Coons patches can be efficiently detected and split into nondegen- 

477 erate subpatches. This is required to avoid both patch overlap at the bound- 

478 aries of the underlying curve and rendering artifacts (patchiness, smudges, 

479 or wrinkles) due to normal reversal. 

480 We have illustrated applications of these techniques in the open source 

481 vector graphics programming language ASYMPTOTE, which we believe is the 

482 first software to lift T[t;X into 3D. This represents an important milestone for 

483 publication-quality scientific graphing. 
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484 Appendix A. Bezier Approximation of a Sphere 

485 As previously emphasized, although conic sections (quadrics) may be ac- 

486 curately represented by NURBS surfaces, the language of high-end printers, 

487 PostScript, supports only Bezier curves and surfaces. Although PostScript 

488 is only a 2D language, vector graphics projections of Bezier surfaces are 

489 nevertheless possible using tensor-product patch shading and hidden surface 

490 splitting along approximations to the visible surface horizon. 

491 Here we illustrate that a sphere may be approximated to high graphical 

492 accuracy by a Bezier surface with only 8 patches, one for each octant, fol- 

493 lowing a procedure suggested in Ref. [18]. The patch describing an octant 

494 is degenerate at the pole: two of the nodes are placed there, with the other 

495 two placed along the equator, 90° apart in longitude. 

496 Following Knuth, a unit quarter circle is approximated in ASYMPTOTE 

497 "with less than 0.06% error" [12], using the control points {(1, 0), (1, a), (a, 1), (0, 1)}, 

498 where a = |(v2 — 1). This value of a is determined by requiring that the 

499 third-order Bezier midpoint lie on the unit circle at (l/\/2, l/\/2)- (Other 

500 methods of approximating circular arcs by Bezier curves have been described 

501 inRefs. [5], 0, and [2I].) 

502 The above prescription immediately determines the three circular arcs 

503 describing the patch boundary for a unit spherical octant. Let us place 

504 Poo at (1,0,0), Po3 = Pi3 = P23 = P33 at (0,0,1), and P30 at (0,1,0). 

505 The remaining control points {Pn, P12, P21, P22} are chosen to make the 

506 surface nearly spherical and the interface with adjacent octants smooth (have 

507 continuous first derivatives at the patch boundaries). The point Pn is chosen 

508 (on the tangent plane at a; = 1) to be the vector sum Pio + Poi — Poo = 

509 (l,a, 0) -|- (1,0, a) — (1,0,0) = (l,a, a). We also require that the triangle 

510 in the x-y plane formed by the origin and the projections of P12 onto the 

511 x-y plane and the x axis is similar to the corresponding triangle for Pn. 

512 This implies that P12 = (a,a^,l). Similarly, we determine P22 = (a^,a, 1) 

513 and P21 = (a, l,a). The final Bezier patch and resulting approximation to 

514 a unit sphere, with the control point mesh shown in blue, are illustrated 



515 in Figure A. 12 We found numerically that the radius of this approximate 

516 sphere, generated with a 12 x 7 control point mesh, varies by less than 0.052%, 

517 well below the tolerance 0.1% to which Figure 8 of Ref. [2D] was drawn using 

518 a much finer 22 x 13 control point mesh. 
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(1,0, a). 



(1,0,0) 




(l,a,0) 




Figure A. 12: Bezier approximation to a unit sphere. The red dots indicate control points. 
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